#Code built for Gabriel Fulton thesis in order to import raw text files
#Imported from text for processing, graphing, and analysis
#packages: quantmod, DT (for seperate tables)
#packages?: TTR, xts, zoo

#Set wd to whichever folder currently has the files which will be analyzed
#Surface
setwd("C:/Users/gfult/Dropbox/Salt/InProgress/Spectro Analysis/08Jun2018/BeetSpectraRedux/PSR+3500_1596061/2018_Jun_08")
#Home PC
#setwd("C:/Users/Gabriel Fulton/Dropbox/Salt/InProgress/Spectro Analysis/08Jun2018/BeetSpectraRedux/PSR+3500_1596061/2018_Jun_08")

#Read from .sed file based on spacing using fwf, "View(anysample)" for check
S1=read.fwf("Jun8beet_00061.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S2=read.fwf("Jun8beet_00062.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S3=read.fwf("Jun8beet_00063.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S4=read.fwf("Jun8beet_00064.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S5=read.fwf("Jun8beet_00065.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S6=read.fwf("Jun8beet_00066.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S7=read.fwf("Jun8beet_00067.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S8=read.fwf("Jun8beet_00068.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S9=read.fwf("Jun8beet_00069.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S10=read.fwf("Jun8beet_00070.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)

#Graph all individual readings (Currently can't fit 10...) [layout function cool af]
#x11() makes a serpate window for graph
#x11()
#par(mfrow=c(2,5))
#plot(S1)
#plot(S2)
#plot(S3)
#plot(S4)
#plot(S5)
#plot(S6)
#plot(S7)
#plot(S8)
#plot(S9)
#plot(S10)


#Average these spectral readings with each other
beetdish100=(S1+S2+S3+S4+S5+S6+S7+S8+S9+S10)/10

#Plot the cumulative average
x11()
par(mfrow=c(1,1))
plot(beetdish100, main="")

#Divided by highest value such that max is now 1
reflectance=beetdish100$Reflectance/max(beetdish100$Reflectance)
#x11()
#par(mfrow=c(1,1))

beetdish100DividedByMax=cbind.data.frame(350:2500,reflectance)
#plot(beetdishsalt5DividedByMax)

#Correlation
#Truncating data for correlation, 350-1400
#Analysis of key points (Salt) - 1490, 1990
#this covers 350 to 1400 nm wavelength
tbeetdish100=beetdish100[1:1051,1:2]
cor(tsalt$Reflectance,tbeetdish100$Reflectance)

#Divide out the water and dish
cor(tsalt$Reflectance,tbeetdish100$Reflectance/twaterdish$Reflectance)

#Taking natural log
#maybe divide out beet juice equivalence instead of waterdish0
cor(log(tsalt$Reflectance),log(tbeetdish100$Reflectance/twaterdish$Reflectance))

#Taking natural log as step 1
#maybe divide out beet juice equivalence instead of waterdish0
cor(log(tsalt$Reflectance),log(tbeetdish100$Reflectance)/log(twaterdish$Reflectance))

#Grab desired POI y values (at wavelength = 986,1001,1030,1058,1193,1256,1457,1632)
POIbeetdish100=rbind(beetdish100[637,],beetdish100[652,],beetdish100[681,],beetdish100[709,],beetdish100[844,],beetdish100[907,],beetdish100[1108,],beetdish100[1283,])
POIbeetdish100

#Salt derived POI
POI2beetdish100=rbind(beetdish100[129,],beetdish100[509,],beetdish100[557,],beetdish100[637,],beetdish100[744,],beetdish100[874,],beetdish100[921,],beetdish100[1308,])
POI2beetdish100

#Prints to notepad to copy to excel
write.table(beetdish100,"C:/Users/gfult/Desktop/Upload/beetdish100.txt",sep="\t")

#EVERYTHING BELOW HERE IS WORK IN PROGRESS FOR AUTOMATION OF PEAK FINDING

#Apparently embedded findPeaks is no good so adding this function:
find_peaks <- function (x, m = 3){
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
    z <- i - m + 1
    z <- ifelse(z > 0, z, 1)
    w <- i + m + 1
    w <- ifelse(w < length(x), w, length(x))
    if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
  })
  pks <- unlist(pks)
  pks
}
#Where m is number of points on either side of the peak
#the function can also be used to find local minima of any sequential vector x via find_peaks(-x)

#THIS AREA IS MESSED UP AND NEEDS WORK
#Peak and valley analysis (NOT SURE OF WHAT m IS APPROPRIATE, originally 3)
peaks=find_peaks(Stotal$Reflectance, m=30)
valleys=find_peaks(-Stotal$Reflectance, m=30)

#peaks=findPeaks(Stotal$Reflectance,thresh=0.05)
#valleys=findValleys(Stotal$Reflectance,thresh=0.05)

#Finds corresponding reflectance values for found wavelength peaks and valleys
ppeaks=Stotal[Stotal$Wavelength%in%peaks,]
pvalleys=Stotal[Stotal$Wavelength%in%valleys,]


#Highlight of said peak and valleys with corresponding wavelengths/severity
points(ppeaks,col="red",pch=19)
points(pvalleys,col="blue",pch=19)

#Creation of tables of corresponding data:
#Peaks:
View(ppeaks, "Peaks")
#Valleys:
View(pvalleys, "Valleys")

#For Table writing:
#write.xlsx(x, file, sheetName="Sheet1")

